#include "phys.h"

/* ---------------------------------------------------------------- */
/* ---------------------------------------------------------------- */

Phys_var::Phys_var() {

  dx	=  0;
  dy   	=  0;
  dt 	=  0;
  Jn 	=  0;
  Jp 	=  0;
  D	=  0;
  alloy_x = 0;
  alloy_y = 0;
  phi	=  0;
  phi_n =  0;
  phi_p =  0;
  qphin =  0;
  qphip =  0;
  a =  0;
  c =  0;
  doping = 0;
  elem 	=  0;
  qptr_e =  0;
  qptr_h =  0;

}

/* ---------------------------------------------------------------- */
/* ---------------------------------------------------------------- */

Phys_var::Phys_var(const Phys_var & op) {

  dx	=  op.dx;
  dy   	=  op.dy;
  dt 	=  op.dt;
  Jn 	=  op.Jn;
  Jp 	=  op.Jp;
  D	=  op.D;
  alloy_x = op.alloy_x;
  alloy_y = op.alloy_y;
  phi	=  op.phi;
  phi_n =  op.phi_n;
  phi_p =  op.phi_p;
  qphin =  op.qphin;;
  qphip =  op.qphip;;
  a =  op.a;
  c =  op.c;
  doping = op.doping;
  elem 	=  op.elem;
  if(op.qptr_e) {
    qptr_e = new Qntm_var;
    *qptr_e = *(op.qptr_e);
  }
  else qptr_e = 0;
  if(op.qptr_h) {
    qptr_h = new Qntm_var;
    *qptr_h = *(op.qptr_h);
  }
  else qptr_h = 0;

}

/* ---------------------------------------------------------------- */
/* ---------------------------------------------------------------- */

Phys_var Phys_var::operator=(const Phys_var & op) {

  dx	=  op.dx;
  dy   	=  op.dy;
  dt 	=  op.dt;
  Jn 	=  op.Jn;
  Jp 	=  op.Jp;
  D	=  op.D;
  alloy_x = op.alloy_x;
  alloy_y = op.alloy_y;
  phi	=  op.phi;
  phi_n =  op.phi_n;
  phi_p =  op.phi_p;
  qphin =  op.qphin;;
  qphip =  op.qphip;;
  a =  op.a;
  c =  op.c;
  doping = op.doping;
  elem 	=  op.elem;
  if(op.qptr_e) {
    qptr_e = new Qntm_var;
    *qptr_e = *(op.qptr_e);
  }
  else qptr_e = 0;
  if(op.qptr_h) {
    qptr_h = new Qntm_var;
    *qptr_h = *(op.qptr_h);
  }
  else qptr_h = 0;

  return * this;

}
/* ---------------------------------------------------------------- */
/* ---------------------------------------------------------------- */
	
void Phys_var::set ( double dx_inp,double dy_inp,
		     double dt_inp,double alloy_x_inp,double alloy_y_inp,
		     double doping_inp,double bias,
		     Alloy_var * ref,Qntm_var * ptr_e,Qntm_var * ptr_h ) {

  double Eg,chi,Vbi,p,n,Nc,Nv,y,E_ion,up,down,traps;


  dx	=  dx_inp;
  dy	=  dy_inp;
  dt 	=  dt_inp;
  Jn 	=  0;
  Jp 	=  0;
  D	=  0;
  a     =  0;
  c     =  0;
  doping = doping_inp;
  alloy_x = alloy_x_inp;
  alloy_y = alloy_y_inp;
  elem	=  ref;

  Eg 	= elem->Eg(alloy_x,alloy_y);
  chi 	= elem->chi(alloy_x,alloy_y);
  Nc 	= elem->Nc(alloy_x,alloy_y);
  Nv 	= elem->Nv(alloy_x,alloy_y);
  E_ion = elem->Eion(doping);
		
  down	= -Eg;
  up   	= Eg;

  do {	/* Vbi = -phi-chi+phi_n = Ec-Ef */
    Vbi	= (up+down)/2.0;
    n = Nc*FD(-Vbi/VT);
    p = Nv*FD((Vbi-Eg)/VT);
    if(elem->Ntrap()>0 || elem->Ptrap()>0 )
      traps = elem->Ntrap()*GX((Vbi-elem->Etrap_e())/VT,1.0) -elem->Ptrap()*GX((Vbi-Eg-elem->Etrap_h())/VT,1.0);
    else traps = 0.0;
    if(doping>=0) y = p - n + doping*GD(-(Vbi-E_ion)/VT) - traps;
    else 	  y = p - n + doping*GA(-(Vbi-Eg+E_ion)/VT) - traps;
    if(y<0) down = Vbi;
    else    up = Vbi;
  } while( (up-down>1e-15) );

  phi   = bias + ( -Vbi + elem->chi(0,0) - chi);
  phi_n = phi  + chi + Vbi;
  phi_p = phi_n;
  qphin = phi_n;
  qphip = phi_n;

  if(ptr_e!=NULL) {
    qptr_e = new Qntm_var;
    *qptr_e = *ptr_e;
  }
  else  qptr_e = 0;

  if(ptr_h!=NULL) {
    qptr_h = new Qntm_var;
    *qptr_h = *ptr_h;
  }
  else qptr_h = 0;

}

/* ---------------------------------------------------------------- */
/* ---------------------------------------------------------------- */
	
void Phys_var::set_potential ( double bias ) {

  double Eg,chi,Vbi,p,n,Nc,Nv,y,E_ion,up,down,traps;

  Eg 	= elem->Eg(alloy_x,alloy_y);
  chi 	= elem->chi(alloy_x,alloy_y);
  Nc 	= elem->Nc(alloy_x,alloy_y);
  Nv 	= elem->Nv(alloy_x,alloy_y);
  E_ion = elem->Eion(doping);
		
  down	= -Eg;
  up   	= Eg;

  do {		/* Vbi = -phi-chi+phi_n = Ec-Ef */
    Vbi	= (up+down)/2.0;
    n = Nc*FD(-Vbi/VT);
    p = Nv*FD((Vbi-Eg)/VT);
    if(elem->Ntrap()>0 || elem->Ptrap()>0 )
      traps = elem->Ntrap()*GX((Vbi-elem->Etrap_e())/VT,1.0) -elem->Ptrap()*GX((Vbi-Eg-elem->Etrap_h())/VT,1.0);
    else traps = 0.0;
    if(doping>=0) y = p - n + doping*GD(-(Vbi-E_ion)/VT) -traps;
    else 	  y = p - n + doping*GA(-(Vbi-Eg+E_ion)/VT)-traps;
    if(y<0) down = Vbi;
    else    up = Vbi;
  } while( (up-down>1e-15) );

  phi   = bias + ( -Vbi + elem->chi(0,0) - chi);
  phi_n = phi  + chi + Vbi;
  phi_p = phi_n;
  qphin = phi_n;
  qphip = phi_n;

}

/* ---------------------------------------------------------------- */
/* ---------------------------------------------------------------- */


void Phys_var::show_optical_bandgap(FILE * fp) {
  int j;
  const double refractive_index = 1.0;
  for(j=0;j<5;j++) {
    if(qptr_e!=NULL && j<qptr_e->ret_num() && qptr_h!=NULL && j<qptr_h->ret_num())
      fprintf(fp,"%le  %le  ",qptr_e->ret_eval(j).energy+qptr_h->ret_eval(j).energy,
	      refractive_index*1241.54/(qptr_e->ret_eval(j).energy+qptr_h->ret_eval(j).energy));
    else fprintf(fp,"0.0e0  ");
  }
}

void Phys_var::show_lev_e(FILE * fp) {
  int j;
  for(j=0;j<10;j++) {
    if(qptr_e!=NULL && j<qptr_e->ret_num()) fprintf(fp,"%le  ",qptr_e->ret_eval(j).energy);
    else fprintf(fp,"0.0e0  ");
  }
}

void Phys_var::show_lev_h(FILE * fp) {
  int j;
  for(j=0;j<10;j++) {
    if(qptr_h!=NULL && j<qptr_h->ret_num()) fprintf(fp,"%le  ",-qptr_h->ret_eval(j).energy);
    else fprintf(fp,"0.0e0  ");
  }
}

void Phys_var::show_well_e(FILE * fp) {
  int j;
  for(j=0;j<10;j++) {
    if(qptr_e!=NULL && j<qptr_e->ret_num()) fprintf(fp,"%le  ",qptr_e->prob[j]);
    else fprintf(fp,"0.0e0  ");
  }
}

void Phys_var::show_well_h(FILE * fp) {
  int j;
  for(j=0;j<10;j++) {
    if(qptr_h!=NULL && qptr_h->ret_num()>j) fprintf(fp,"%le  ",qptr_h->prob[j]);
    else fprintf(fp,"0.0e0  ");
  }
}

  
/* ---------------------------------------------------------------- */
/* ---------------------------------------------------------------- */

 
double Phys_var::n_sum(double Ein, double Efin, double Ec, double Ef) {

  double sum=0,h,E,E0;

  E0 = (Ein>=Ec)? Ec : Ein;
  if(Efin<E0) nrerror("N_SUM: Efin<Ein");
  h = (Efin-E0)/100.0;
  E = Ein;
  while(E<=Efin) {
    sum += n(E,Ec,Ef)*h;
    E+=h;
  }

  return sum;

}

/* ---------------------------------------------------------------- */
/* ---------------------------------------------------------------- */

void Phys_var::set_c(double a_inp, double x_inp, double y_inp) {

  double a0,c0,poisson_ratio;

  if(a_inp > 0) {
    a0 = elem->a(x_inp,y_inp);
    c0 = elem->c(x_inp,y_inp);
    poisson_ratio = elem->poisson_ratio(x_inp,y_inp);  
    c = c0*(1.0-2.0*(poisson_ratio/(1.0-poisson_ratio))*((a_inp-a0)/a0));
  }
  else c = 0.0;

}

/* ---------------------------------------------------------------- */
/* ---------------------------------------------------------------- */

double Phys_var::Opt_Rec() {

  double ovrlp;
  double T_optical = 1e-8;
  double rec_base = ((n()>p())?p():n())/T_optical;

  if(qptr_e!=0 && qptr_h!=0 && qptr_e->prob && qptr_h->prob ) {
    ovrlp = qptr_e->prob[0]*qptr_h->prob[0];
    rec_base = rec_base*(1.0+1000.0*ovrlp );
  }
  return rec_base;

}

/* ---------------------------------------------------------------- */
/* ---------------------------------------------------------------- */

double Phys_var::QW_Rec() {

  double T_capture = 1e-9;
  double T_emission = 1e-7;
  double rec_base;

  rec_base = (n()/T_capture)-nqw()/T_emission;
  return rec_base;

}

// ---------------------------------------------------------------
// ---------------------------------------------------------------

// Calculate the density of electrons in the quantum well
// It returns a very low value if the quantum well is not present.

double Phys_var::nqw(){ 
	int i;
    double tmp = 0;
    if(!qptr_e) return 1e-40;
    for(i=0;i<qptr_e->ret_num();i++) tmp += me()*N2DEG*log(1.0+exp((-qphin-qptr_e->ret_eval(i).energy)/VT))*qptr_e->prob[i];
    if(tmp<1e-40) tmp = 1e-40;
    return tmp;
  }

// ---------------------------------------------------------------
// ---------------------------------------------------------------

// Calculate the density of holes in the quantum well

double Phys_var::pqw() { 
    int i;
    double tmp = 0;
    if(!qptr_h) return 1e-40;
    for(i=0;i<qptr_h->ret_num();i++) tmp += mh()*N2DEG*log(1.0+exp((qptr_h->ret_eval(i).energy+qphip)/VT))*qptr_h->prob[i];
    if(tmp<1e-40) tmp = 1e-40;
    return tmp; 
  }